{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "遥感影像的裁剪和拼接\n",
    "  注意坐标原点在左上角，x轴正方向向右，y轴正方向向下"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "from osgeo import gdal\n",
    "import os\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "extent: 102.15087890625 103.68896484375 24.36767578125 26.56494140625\n"
     ]
    }
   ],
   "source": [
    "path_img=\"D:\\Lenovo\\Desktop\\云南大学\\开源地理信息系统\\output\\昆明市_高程_Level_14.tif\"\n",
    "dset=gdal.Open(path_img)\n",
    "geo_trans=dset.GetGeoTransform()\n",
    "geo_trans\n",
    "x_min,y_max=geo_trans[0],geo_trans[3] #x_min,y_max分别对应第一个和第四个索引\n",
    "x_max=x_min+geo_trans[1]*dset.RasterXSize\n",
    "y_min=y_max+geo_trans[5]*dset.RasterYSize #也可写为geo_trans[-1]*dset.RasterYSize，-1代表最后一位\n",
    "print('extent:',x_min,x_max,y_min,y_max)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_min_sub,x_max_sub=102.5,103\n",
    "y_min_sub,y_max_sub=24.5,25.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "775.4666666666667 254.22222222222223\n"
     ]
    }
   ],
   "source": [
    "#求给定范围左上角坐标的图像坐标\n",
    "\n",
    "#计算X方向上（列方向）的起点坐标\n",
    "col_start_subs=(x_min_sub-x_min)/geo_trans[1]\n",
    "#计算Y方向上（行方向）的起点坐标\n",
    "row_start_subs=(y_max_sub-y_max)/geo_trans[5]\n",
    "print(row_start_subs,col_start_subs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "775.4666666666667 254.22222222222223\n",
      "102.49969482421875 25.500640869140625\n"
     ]
    }
   ],
   "source": [
    "#对裁剪影像上的起点坐标进行整数化\n",
    "col_start_subs_update,row_start_subs_update=round(col_start_subs),round(row_start_subs)\n",
    "print(row_start_subs,col_start_subs)\n",
    "#使图像上的坐标由实际地理坐标转为像元中点值坐标\n",
    "x_min_subs_update=geo_trans[0]+geo_trans[1]*col_start_subs_update\n",
    "y_max_subs_update=geo_trans[3]+geo_trans[5]*row_start_subs_update\n",
    "print(x_min_subs_update,y_max_subs_update)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "裁剪范围在图像上的图像尺寸"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "364.31111111111113 728.6444444444444\n",
      "364 729\n"
     ]
    }
   ],
   "source": [
    "x_size_subs=(x_max_sub-x_min_subs_update)/geo_trans[1]\n",
    "y_size_subs=(y_min_sub-y_max_subs_update)/geo_trans[5]\n",
    "print(x_size_subs,y_size_subs)\n",
    "x_size_subs_update,y_size_subs_update=round(x_size_subs),round(y_size_subs)\n",
    "print(x_size_subs_update,y_size_subs_update)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "102.99957275390625 24.49951171875\n",
      "更新后的裁剪范围 [102.49969482421875, 102.99957275390625, 24.49951171875, 25.500640869140625]\n"
     ]
    }
   ],
   "source": [
    "x_max_subs_update=x_min_subs_update+geo_trans[1]*x_size_subs_update\n",
    "y_min_subs_update=y_max_subs_update+geo_trans[5]*y_size_subs_update\n",
    "print(x_max_subs_update,y_min_subs_update)\n",
    "print('更新后的裁剪范围',[x_min_subs_update,x_max_subs_update,y_min_subs_update,y_max_subs_update])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1600, 1120)"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#裁剪影像的地理转换参数\n",
    "geotrans_subs=[x_min_subs_update,20,0,y_max_subs_update,0,-20]\n",
    "#输出数组\n",
    "img_array=dset.ReadAsArray()\n",
    "img_array.shape\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(729, 364)"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "img_array_subs=img_array[\n",
    "                         row_start_subs_update:row_start_subs_update+y_size_subs_update,\n",
    "                         col_start_subs_update:col_start_subs_update+x_size_subs_update]\n",
    "img_array_subs.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "写出裁剪影像"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "expected array of dim 2",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[1;32md:\\Lenovo\\Desktop\\云南大学\\开源地理信息系统\\Code\\Chapter5\\20231016.ipynb 单元格 13\u001b[0m line \u001b[0;36m1\n\u001b[0;32m      <a href='vscode-notebook-cell:/d%3A/Lenovo/Desktop/%E4%BA%91%E5%8D%97%E5%A4%A7%E5%AD%A6/%E5%BC%80%E6%BA%90%E5%9C%B0%E7%90%86%E4%BF%A1%E6%81%AF%E7%B3%BB%E7%BB%9F/Code/Chapter5/20231016.ipynb#X15sZmlsZQ%3D%3D?line=8'>9</a>\u001b[0m \u001b[39m# for i in range(dset.RasterCount):\u001b[39;00m\n\u001b[0;32m     <a href='vscode-notebook-cell:/d%3A/Lenovo/Desktop/%E4%BA%91%E5%8D%97%E5%A4%A7%E5%AD%A6/%E5%BC%80%E6%BA%90%E5%9C%B0%E7%90%86%E4%BF%A1%E6%81%AF%E7%B3%BB%E7%BB%9F/Code/Chapter5/20231016.ipynb#X15sZmlsZQ%3D%3D?line=9'>10</a>\u001b[0m outband\u001b[39m=\u001b[39mdset\u001b[39m.\u001b[39mGetRasterBand(\u001b[39m1\u001b[39m)\n\u001b[1;32m---> <a href='vscode-notebook-cell:/d%3A/Lenovo/Desktop/%E4%BA%91%E5%8D%97%E5%A4%A7%E5%AD%A6/%E5%BC%80%E6%BA%90%E5%9C%B0%E7%90%86%E4%BF%A1%E6%81%AF%E7%B3%BB%E7%BB%9F/Code/Chapter5/20231016.ipynb#X15sZmlsZQ%3D%3D?line=10'>11</a>\u001b[0m outband\u001b[39m.\u001b[39;49mWriteArray(img_array_subs[\u001b[39m0\u001b[39;49m])\n\u001b[0;32m     <a href='vscode-notebook-cell:/d%3A/Lenovo/Desktop/%E4%BA%91%E5%8D%97%E5%A4%A7%E5%AD%A6/%E5%BC%80%E6%BA%90%E5%9C%B0%E7%90%86%E4%BF%A1%E6%81%AF%E7%B3%BB%E7%BB%9F/Code/Chapter5/20231016.ipynb#X15sZmlsZQ%3D%3D?line=11'>12</a>\u001b[0m dset_subs\u001b[39m=\u001b[39m\u001b[39mNone\u001b[39;00m\n",
      "File \u001b[1;32mc:\\Users\\Lenovo\\.conda\\envs\\myenv\\lib\\site-packages\\osgeo\\gdal.py:4740\u001b[0m, in \u001b[0;36mBand.WriteArray\u001b[1;34m(self, array, xoff, yoff, resample_alg, callback, callback_data)\u001b[0m\n\u001b[0;32m   4734\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mWriteArray\u001b[39m(\u001b[39mself\u001b[39m, array, xoff\u001b[39m=\u001b[39m\u001b[39m0\u001b[39m, yoff\u001b[39m=\u001b[39m\u001b[39m0\u001b[39m,\n\u001b[0;32m   4735\u001b[0m                resample_alg\u001b[39m=\u001b[39mgdalconst\u001b[39m.\u001b[39mGRIORA_NearestNeighbour,\n\u001b[0;32m   4736\u001b[0m                callback\u001b[39m=\u001b[39m\u001b[39mNone\u001b[39;00m,\n\u001b[0;32m   4737\u001b[0m                callback_data\u001b[39m=\u001b[39m\u001b[39mNone\u001b[39;00m):\n\u001b[0;32m   4738\u001b[0m     \u001b[39mfrom\u001b[39;00m \u001b[39mosgeo\u001b[39;00m \u001b[39mimport\u001b[39;00m gdal_array\n\u001b[1;32m-> 4740\u001b[0m     \u001b[39mreturn\u001b[39;00m gdal_array\u001b[39m.\u001b[39;49mBandWriteArray(\u001b[39mself\u001b[39;49m, array, xoff, yoff,\n\u001b[0;32m   4741\u001b[0m                                       resample_alg\u001b[39m=\u001b[39;49mresample_alg,\n\u001b[0;32m   4742\u001b[0m                                       callback\u001b[39m=\u001b[39;49mcallback,\n\u001b[0;32m   4743\u001b[0m                                       callback_data\u001b[39m=\u001b[39;49mcallback_data)\n",
      "File \u001b[1;32mc:\\Users\\Lenovo\\.conda\\envs\\myenv\\lib\\site-packages\\osgeo\\gdal_array.py:465\u001b[0m, in \u001b[0;36mBandWriteArray\u001b[1;34m(band, array, xoff, yoff, resample_alg, callback, callback_data)\u001b[0m\n\u001b[0;32m    461\u001b[0m \u001b[39m\u001b[39m\u001b[39m\"\"\"Pure python implementation of writing a chunk of a GDAL file\u001b[39;00m\n\u001b[0;32m    462\u001b[0m \u001b[39mfrom a numpy array.  Used by the gdal.Band.WriteArray method.\"\"\"\u001b[39;00m\n\u001b[0;32m    464\u001b[0m \u001b[39mif\u001b[39;00m array \u001b[39mis\u001b[39;00m \u001b[39mNone\u001b[39;00m \u001b[39mor\u001b[39;00m \u001b[39mlen\u001b[39m(array\u001b[39m.\u001b[39mshape) \u001b[39m!=\u001b[39m \u001b[39m2\u001b[39m:\n\u001b[1;32m--> 465\u001b[0m     \u001b[39mraise\u001b[39;00m \u001b[39mValueError\u001b[39;00m(\u001b[39m\"\u001b[39m\u001b[39mexpected array of dim 2\u001b[39m\u001b[39m\"\u001b[39m)\n\u001b[0;32m    467\u001b[0m xsize \u001b[39m=\u001b[39m array\u001b[39m.\u001b[39mshape[\u001b[39m1\u001b[39m]\n\u001b[0;32m    468\u001b[0m ysize \u001b[39m=\u001b[39m array\u001b[39m.\u001b[39mshape[\u001b[39m0\u001b[39m]\n",
      "\u001b[1;31mValueError\u001b[0m: expected array of dim 2"
     ]
    }
   ],
   "source": [
    "driver=gdal.GetDriverByName('GTiff')\n",
    "path_img_out=\"D:\\Lenovo\\Desktop\\云南大学\\开源地理信息系统\\output\\subset_20231018.tif\"\n",
    "#创建数组驱动\n",
    "dset_subs=driver.Create(path_img_out,xsize=x_size_subs_update,ysize=y_size_subs_update,bands=dset.RasterCount,eType=gdal.GDT_Int16)\n",
    "\n",
    "#定义地理转换\n",
    "dset_subs.SetGeoTransform(geotrans_subs)\n",
    "dset_subs.SetProjection(dset.GetProjection())\n",
    "for i in range(dset.RasterCount):\n",
    "    outband=dset.GetRasterBand(i+1)\n",
    "    outband.WriteArray(img_array_subs[i])\n",
    "dset_subs=None\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
